# TODO: Add comment
# 
# Author: broken
###############################################################################
library(raster)
library(base)
library(gstat)

b=as(extent(102.144996643067,109.469429016114,8.56333160400391,23.3927307128908),"SpatialPolygons")

rasterFile=raster("C:/Users/Forever/Desktop/dem_10km.tif")
rasterFile=crop(rasterFile,b)
writeRaster(rasterFile,filename="C:/Users/Forever/Desktop/dem10.tif",format="GTiff")


dem=raster("C:/Users/Forever/Desktop/dem.tif")
aot=raster("C:/Users/Forever/Desktop/aot.tif")
lc=raster("C:/Users/Forever/Desktop/lc.tif")
nl=raster("C:/Users/Forever/Desktop/nl.tif")



olddem=raster("C:/Users/Forever/Desktop/olddem.tif")
newdem=raster("C:/Users/Forever/Desktop/newdem.tif")

oldvalue=values(olddem)
newvalue=values(newdem)

correl=cor(oldvalue,newvalue)
rmse=sqrt(sum((oldvalue-newvalue)^2)/ncell(oldvalue))
re=sum(abs(oldvalue-newvalue)/oldvalue)/ncell(oldvalue)
